#include "StVenKir.hpp"

// variational formulation language
#include <feel/feelvf/vf.hpp>

#include <svk_cl.cpp>
#include <research/life-workspace/Fluid-Structure/functionSup.cpp>

namespace Feel
{


    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    void
    StVenantKirchhoff<Dim, Order, Entity>::updateResidual( const vector_ptrtype& X, vector_ptrtype& R )
    {
#if 0
        using namespace Feel::vf;

        if (M_LookComment) std::cout << "[StVenantKirchhoff] : updateResidual start\n";

        this->changeRepository();

        mesh_ptrtype mesh = M_Xh->mesh();

        element_struct_type u( M_Xh, "u" ); u = *X;
        element_struct_type v( M_Xh, "v" );
        element_struct_type force( M_Xh, "f" );


        //    AUTO( eta, 1000*Px()*( Px() -4 ) /**(Px()-2.5)*/  ) ;
        //AUTO( eta, 1000*(Px()-0.4)*( Px() -0.45 )*chi( Px()>=0.4 && Px()<=0.45 ) ) ;
        //AUTO( eta, 0.2*(Px()-0.4)*( Px() -0.45 )*chi( Px()>=0.4 && Px()<=0.45 ) ) ;
        //AUTO( eta, 0.2*(Px()-1.4)*( Px() -1.45 )*chi( Px()>=1.4 && Px()<=1.45 ) ) ;
        auto eta = 0.3*sin(M_PI*Px()/4.);

        //auto eta = -0.7*sin(M_bdf_displ_struct->time()*M_PI/(0.7*2.) )*0.3*vf::sin( M_PI*vf::Px()/4.);
        //auto grad_eta = mat<2,2>(  cst(0.), cst(0.),-0.7*sin(M_bdf_displ_struct->time()*M_PI/(0.07*2.) )*0.3*(M_PI/4.)*vf::cos( M_PI*vf::Px()/4.), cst(0.) );

        //force = vf::project(M_Xh, elements(mesh), eta*oneY());

        //Matrice Identitee
        auto Id = oneX()*trans(oneX()) + oneY()*trans(oneY());

        auto Fv = Id + gradv(u);
        auto F = Id + grad(u);
        //auto Fvf = Id + grad_eta; //AUTO( Fvf , Id + gradv(force) );

        //AUTO( Evf, 0.5*(gradv(force)+trans(gradv(force)) ) + 0.5*trans(gradv(force))*gradv(force) );
        //auto Evf = 0.5*(grad_eta+trans(grad_eta) ) + 0.5*trans(grad_eta)*grad_eta;
        auto Ev = 0.5*(gradv(u)+trans(gradv(u)) ) + 0.5*trans(gradv(u))*gradv(u);
        auto E = 0.5*(grad(u)+trans(grad(u)) ) + 0.5*trans(grad(u))*grad(u);

        //auto Svf = lambda*trace(Evf)*Id + 2*mu*Evf;
        auto Sv = lambda*trace(Ev)*Id + 2*mu*Ev;
        auto S = lambda*trace(E)*Id + 2*mu*E;

        //auto dFf = grad_eta; //AUTO( dFf, gradv(force));
        auto dF = gradv(u);
        auto dFF = grad(u);

        //auto dEf = 0.5*(trans(Fvf)*dFf + trans(dFf)*Fvf);
        auto dE = 0.5*(trans(Fv)*dF + trans(dF)*Fv);
        auto dEE = 0.5*(trans(Fv)*dFF + trans(dFF)*Fv);

        //auto dSf = lambda*trace(dEf)*Id + 2*mu*dEf;
        auto dS = lambda*trace(dE)*Id + 2*mu*dE;
        auto dSS = lambda*trace(dEE)*Id + 2*mu*dEE;
        /*
        //force volumique :
        AUTO(f,vec(cst(0.),cst(0.)) );
        */
        //--------------------------------------------------------------------------------------------------//

        form1( M_Xh, R, _init=true ) =
            integrate( elements(mesh), _Q<4*(Order)>(),
                       trace((Fv*Sv)*trans(grad(v))) );

        double penalisation_bc=this->application()->vm()["bccoeff"].template as<double>();
#if 0
        form1( M_Xh, R ) +=
            //integrate( elements(mesh), _Q<2*Order>(),
            integrate( markedfaces(mesh,mesh->markerName("Force")), _Q<2*Order>(),
                       - trans(eta*oneY())*id(v) );
#endif


#if 1

        auto clDef = CL_DEF_MECASOL(M_time);

        if (M_weakCL)
            {


                ForEachCL( clDef,cl::dirichlet_vec,
                           form1( M_Xh, R ) +=
                           integrate( markedfaces(mesh,mesh->markerName(PhysicalName)), _Q<4*Order>(),
                                      //terme classique
                                      - trans(Fv*Sv*N())*id(v)
                                      //terme de symetrisation de la jacobienne
                                      - trans((dFF*Sv + Fv*dSS)*N())*idv(u)
                                      + penalisation_bc*trans(idv(u))*id(v)/hFace()
                                      + trans((dFF*Sv + Fv*dSS)*N())*Expression
                                      - penalisation_bc*trans( Expression )*id(v)/hFace()
                                      )
                           );
            }
#endif

        //il manque les neumann ( vide car cond libre)

#if 0
        form1( M_Xh, R ) +=
            integrate( markedfaces(mesh,mesh->markerName("Fixe")), _Q<4*Order>(),
                       //terme classique
                       - trans(Fv*Sv*N())*id(v)
                       //terme de symetrisation de la jacobienne
                       - trans((dFF*Sv + Fv*dSS)*N())*idv(u)
                       //- trans((dF*Sv + Fv*dS)*N())*id(v)
                       // - trans(F*S*N())*idv(u)
                       + penalisation_bc*trans(idv(u))*id(v)/hFace() );
#endif

#if 1
#if 0
        form1( M_Xh, R ) +=
            integrate( markedfaces(mesh,mesh->markerName("Force")), _Q<4*Order>(),
                       //terme classique
                       - trans(Fv*Sv*N())*id(v)
                       //terme de symetrisation
                       - trans((dFF*Sv + Fv*dSS)*N())*idv(u)
                       //- trans((dF*Sv + Fv*dS)*N())*id(v)
                       //- trans(F*S*N())*idv(u)
                       + penalisation_bc*trans(idv(u))*id(v)/hFace() );

        form1( M_Xh, R ) +=
            integrate( markedfaces(mesh,mesh->markerName("Force")), _Q<4*Order>(),
                       //terme de symetrisation
                       + trans((dFF*Sv + Fv*dSS)*N())*idv(force)
                       //+ trans((dFf*Svf + Fvf*dSf)*N())*id(v)
                       // - trans(F*S*N())*idv(u)
                       //!!!- penalisation_bc*trans(idv(force)/*idv(u)*/)*id(v)/hFace() );
                       - penalisation_bc*trans( idv(force)/*eta*oneY()*/)*id(v)/hFace() );
#endif
#else
        form1( M_Xh, R ) +=
            integrate( markedfaces(mesh,mesh->markerName("Force")), _Q<4*Order>(),
                       + trans(Fvf*Svf*N())*id(v) );

#endif


#if 1
        double Newmark_gamma=0.5;
        double Newmark_beta=1/4.;

        auto cst_acc_displ = 1./(Newmark_beta*std::pow(M_bdf_displ_struct->timeStep(),2));
        auto cst_acc_vel = 1./(Newmark_beta*M_bdf_displ_struct->timeStep());
        auto cst_acc_acc = ((1./(2*Newmark_beta))-1);

        element_struct_type buzz2( M_Xh, "u" );
        buzz2.zero();
        buzz2.add(cst_acc_displ,u);
        buzz2.add(-cst_acc_displ, M_bdf_displ_struct->unknown(0));
        buzz2.add(-cst_acc_vel, M_bdf_velocity_struct->unknown(0));
        buzz2.add(-cst_acc_acc, M_bdf_acceleration_struct->unknown(0));

#if 1
        form1( M_Xh, R ) +=
            integrate( elements(mesh), _Q<2*Order>(),
                       trans(idv(buzz2))*id(v) );
#endif
        //buzz = vf::project( M_Xh, elements(mesh), cst_vel_displ*(idv(u)-idv(M_bdf_displ_struct->unknown(0)) ) );

        //buzz->zeros();

        //form2( Xh, Xh, A ) +=
        // integrate( elements(mesh), _Q< 2*Order >(),
        //	 + trans(idt(u))*id(v)*M_bdf_fluid->polyDerivCoefficient(0));
#endif

        R->close();

        if (!M_weakCL)
            {
                auto dirElem = this->getDirichletElement();
                ForEachCL( clDef,cl::dirichlet_vec,
                           modifVec222(*dirElem,R, vf::vec(vf::cst(0.),vf::cst(0.)) , PhysicalName ) );
            }

        if (M_LookComment) std::cout << "[StVenantKirchhoff] : updateResidual finsh\n";
#endif
    }

} // Feel


// instantiation
template class Feel::StVenantKirchhoff<2,2,Feel::Simplex>;
//template class Feel::StVenantKirchhoff<3,2,Feel::Simplex>;
